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ABSTRACT.  ^A  coupled  analytical/numerical  procedure  for  prediction  of 
solute  transport  in  heterogeneous  media  is  described.  The  procedure  con- 
sists of  an  analytic  solution  of  the  hydraulic  equations,  followed  by  a 
numerical  solution  for  solute  transport  using  the  method  of  characteris- 
tics. The  characteristics  are  determined  by  fourth-order  Runge-Kutta  and 
predictor-corrector  algorithms.  Accuracy  of  solute  transport  calculation 
is  enhanced  by  the  fact  that  fluid  velocity  can  be  directly  obtained  at 
priori  undetermined  points  in  the  flow  field. 


The  solute  transport  process  is  considered  to  be  entirely  advective, 
neglecting  the  effects  of  mechanical  dispersion  and  molecular  diffusion. 
Evidence  is  presented  to  demonstrate  that  purely  advective  processes  in 
both  heterogeneous  and  homogeneous  media  can  produce  large  "apparent  dis- 
persion." Such  dispersion  is  shown  to  be  easily  capable  of  overwhelming 
any  reasonable  estimates  of  dispersion  or  diffusion  based  upon  laboratory 
analyses  of  homogeneous  media.  For  groundwater  contamination  problems,  it 
is  concluded  that  precise  definition  of  the  spatial  variability  of  hydrau- 
lic properties  is  crucial  to  the  accurate  determination  of  the  trajectory 
of  contaminated  waters. 


BACKGROUND.  At  the  scale  of  individual  grains,  the  transport  of  a 
conservative  solute  through  a porous  medium  Is  clearly  an  advective 
phenomenon.  Solute  particles  are  wafted  along  by  fluid  as  it  flows  over 
tortuous  routes  in  the  general  direction  of  the  potential  gradient.  Close 
observation  of  the  movement  of  Initially  adjacent  solute  particles  would 
reveal  their  tendency  to  become  separated.  Contributing  to  the  separation 
one  would  observe:  (a)  random  bifurcation  of  pore  channels,  (b)  a large 

range  of  fluid  velocities  across  individual  pores,  and  (c)  differences  of 
fluid  velocity  from  one  pore  to  another.  To  a very  minor  extent,  pore 
scale  advection  is  supplemented  by  molecular  diffusion. 


For  a fluid  of  nonuniform  concentration  which  saturates  a porous 
medium,  the  separation  of  solute  particles  amounts  to  mixing,  resulting  in 
changes  of  local  solute  concentration.  Buyevich  et  al.  (1969)  noted  that 
the  pore  scale  advective  mixing  process  is  very  much  like  the  mixing  re- 
sulting from  ordinary  fluid  turbulence. 

The  usual  material  continuum  approach  to  porous  media  modeling  defines 
solute  transport  in  terms  of  macroscale  mass  fluxes  (Bachmat  and  Bear, 
1964).  Although  the  entire  process  is  fundamentally  advective  at  the  pore 
scale,  the  macroscale  description  of  advection  can  represent  only  an  aver- 
age transport,  strictly  in  the  direction  of  the  potential  gradient.  The 
pore  scale  mechanisms  listed  above  as  (a)  through  (c)  cannot  be  accounted 
for  by  macroscale  advection  alone. 
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Dispersion  is  an  additional  macroscale  solute  flux  whose  purpose  is  to 
account  for  the  pore  scale  mechanisms  which  cause  mixing.  Dispersive  flux 
is  assumed  to  be  proportional  to  concentration  gradient. 

Harleman  et  al.  (1963),  Klotz  and  Moser  (1974),  and  others  have  con- 
ducted laboratory  experiments  on  columns  of  homogeneous  media  to  determine 
the  relative  magnitudes  of  macroscale  advection  and  dispersion.  Their  re- 
sults are  presented  as  a correlation  between  the  magnitude  of  dispersion, 
the  potential  gradient,  and  the  physical  properties  of  fluid  and  media. 

Many  investigators  (c.g.  Pinder,  1973;  Konikow  and  Bredehoeft,  1974) 
have  applied  macroscale  advection-dispersion  models  to  field  problems. 
Calibration  of  these  models  generally  leads  to  the  assumption  of  dispersive 
fluxes  which  are  orders  of  magnitude  greater  than  would  be  expected  on  the 
basis  of  lab  analyses  of  porous  material  samples.  Gelhar  et  al.  (1979) 
reiterate  the  conclusion  that  this  discrepancy  is  related  to  local  hetero- 
geneity of  porous  medium  hydraulic  properties.  The  experimental  results  of 
Skibitski  and  Robinson  (1963)  substantiate  this  by  illustrating  the 
dominant  effect  of  heterogeneity  on  the  transport  of  dye  in  sand  flumes. 

THEORY.  The  aim  of  this  paper  is  to  demonstrate  a coupled  analyti- 
cal/numerical technique  for  predicting  the  transport  of  conservative 
solutes  in  heterogeneous  media.  The  approach  presumes  that  genuine  disper- 
sion is  negligible  compared  to  true  macroscale  advection  when  that  advec- 
tion fully  accounts  for  heterogeneity  and  nonuniform  flow. 

In  order  to  accurately  determine  the  effect  of  heterogeneity  on  macro- 
scale advection,  an  analytic  solution  for  hydraulic  potential  is  obtained. 
Application  of  Darcy’s  law  yields  an  analytic  expression  for  average 
linear  velocity  which  can  be  evaluated  at  priori  unspecified  points  x. 

As  part  of  the  technique,  medium  properties  are  accounted  for  as 
known  (or  interpolated)  explicit  functions  of  x.  Given  an  accurate  de- 
scription of  flow  field,  streamlines  are  calculated  by  applying  the  method 
of  characteristics.  Advection  is  determined  from  the  rates  of  flow  along 
the  streamlines. 

Example  problems  are  used  to  demonstrate  the  fact  that  genuine  disper- 
sion can  be  easily  overwhelmed  by  the  effects  of  heterogeneity  and  non- 
uniform  flow.  In  each  example,  flow  is  assumed  steady  and  horizontal.  The 
analyses  apply  to  confined  aquifers  and  also  to  phreatic  aquifers  where  the 
Dupuit  assumptions  and  linear  approximation  are  valid. 

Advection-dispersion  equation  (numerical  solution).  The  control 
volume  approach  can  be  used  to  derive  the  advection-dispersion  equation 
(Daly,  1979): 

<f>  ~ - V • (DVC)  + v • 7G  «- 4 + q — (1) 

dt  -p  p 

where: 
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<j>  “ effective  porosity,  dimensionless; 

C * mass  fraction  of  the  pore  fluid;  the  ratio  of  the  mass  of  solute 
in  a given  volume  to  the  total  mass  of  fluid  in  that  volume, 
dimensionless; 

D ■ dispersion  coefficient  tensor,  L2/T; 
v ■ specific  dis charge, L/T; 

C “ contaminant  mass  source/sink  strength  representing  the  exchange 
of  mass  between  fluid  and  porous  matrix,  M/L3T; 
p * fluid  density,  M/L3; 

q = recharged  fluid  mass  strength,  M/L3T  0; 

A 

C * mass  fraction  of  recharged  fluid,  dimensionless. 

Considering  (1),  it  is  clear  that  the  effect  of  fluid  withdrawals  is  not 
felt  directly  through  the  ierms  on  the  RHS  of  the  equation.  However,  with- 
drawals do  affect  the  transport  by  modifying  the  flow  field 
(represented  in  (1)  by  v). 

If  dispersion  is  neglected  compared  to  macroscale  advection,  (1)  can 
be  written,  for  the  simple  case  of  no  source/sink  terms,  as: 


3C 

8t 


+ v 


3C  . 3C 
— + v 
3x  y 


ay 


(2) 


Application  of  the  method  of  characteristics  transforms  (2)  into  the  equi- 
valent system  of  ordinary  differential  equations: 


f E f(x,y*t) 

/ 

ft " / = 


(3) 

(4) 


(5) 


Equations  (3)  and  (4)  are  used  to  determine  the  trajectories  (also  called 
the  characteristic  lines)  of  fluid  particles  in  the  flow  field.  Equation 
(5)  is  simply  a statement  of  the  fact  that  in  the  absence  of  sources  of 
sinks  the  concentration  of  fluid  particles  remains  constant.  It  is  impor- 
tant to  note  that  (3),  (4),  and  (5)  are  hot  independent.  Equation  (5)  is 
only  valid  along  the  trajectories  defined  by  the  joint  solution  of  (3)  and 
(4). 


Determination  of  the  trajectories  of  fluid  particles  is  done  by 
numerically  solving  the  linked  system  of  (3)  and  (4).  A fourth  order 
Runge-Kutta  technique  is  used  to  start  the  procedure  which  can  be  continued 
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by  a more  efficient  predictor-corrector  scheme.  The  numerical  solution 
method  starts  at  a point  (x0,y0)  at  time  zero.  The  concentration  is 
defined  at  all  such  points  by  an  initial  condition.  Using  a time  step  At 
successive  points  (xn,  yn)  along  the  trajectory  of  the  particle  which 
began  at  point  (xQ,  y0)  are  obtained.  The  Runge-Kutta  algorithm  for 


accomplishing  this  is: 

Xn+1  “ Xn  + 1 (&1  + 232  + 233  + (6) 

yaH  ” yn  + I (bl  + 2b2  +2b3  + b4>  (7) 

where: 

al  = At  f(xn.yn»tn)  (8) 

bi  - At  g(xn,yn,tn)  (9) 

a2  * At  f(xr  + a3/2,  yn  + bj/2,  tQ  + At/2)  (10) 

b2  * At  g(xQ  + a^/2,  y^  + b3/2,  t^  + At/2)  (11) 

a3  * At  f(xn  + a2/2,  yfl  + b2/2,  tQ  + At/2)  (12) 

bo  « At  g(x  + a2/2,  y + b2/2,  t + At/2)  (13) 

a4  - At  f(xn  + a3,  yn  + b3,  tQ  + At)  (14) 

b4  * At  g(xR  + a3,  yn  + h3,  tft  + At),  (15) 


and  f and  g are  defined  in  Equations  (3)  and  (4). 

Runge-Kutta  algorithms  belong  to  the  set  of  self-starting  numerical 
solution  methods.  Self ^starting  means  that  the  determination  of  all  suc- 
cessive points  (xn,  yn)  requires  only  the  staffing  point  (x0,  yQ). 

In  other  words,  calculation  of  xn  and  yn  depends  only  on  the  known 
values  xn_j  and  ya_j.  The  set  of  non  self-starting  methods  require  the 
values  (xn,  yn)  to  be  given  at  more  than  one  point  along  the  trajec- 
tory. For  example,  a fourth  order  predictor-corrector  algorithm  called 
Milne's  method  requires  the  values  x0,  xlf  x2,  x3,  yG,  y^,  y2,  y3  to 
calculate  successive  values  of  Xq,  and  y^. 
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Beside  the  question  of  starting  values, the  efficiency  of  the  calcula- 
tion procedure  is  an  important  factor  in  selecting  a numerical  method.  It 
turns  out  that  Milne's  algorithm  is  significantly  more  efficient  than  the 
Runge-Kutta  method,  although  both  are  fourth-order  accurate.  One  numerical 
procedure  proposed  in  this  paper  is  that  which  takes  advantage  of  the 
Runge-Kutta  self-starting  feature  and  the  efficiency  of  Milne's  method. 
Given  a starting  point  (xQ,  yQ),  the  Runge-Kutta  procedure  is  used  to 
obtain  (x^,  y^),  (x2,  y2),  (*3,  At  that  stage  the  necessary  starting 

values  are  available  for  Milne's  method  which  is  then  used  to  generate  suc- 
ceeding points. 

Non  self-starting  methods  typically  assume  constant  At,  whereas  self- 
starting methods  allow  for  change  of  At  at  each  time  step.  For  problems  in 
which  the  frequent  change  of  At  is  desirable,  exclusive  use  of  a self- 
starting method,  such  as  the  Runge-Kutta  algorithm,  is  advised. 

Milne's  predictor-corrector  method  consists  of  two  steps.  First,  pre- 
dicted estimates  of  x^  and  yn+i  are  calculated.  Let  these  be  denoted 
x*n+i  and  y*n+i.  Second,  the  predicted  values  are  corrected  to  obtain 
the  final  values  xn+1  and  yn4q  at  the  end  of  a time  step.  The 
algorithm  is  : 


■ vs + i2x;  - 

+ 2xn-2  ^ 

(16) 

»Sfi  ■ Vs + f2y»  - K-i 

+ 2C2]- 

(17) 

then: 

x » x , + 4^  [x*Ii  + ^x'  + 

n+1  n-1  3 1 n+1  n 

(18) 

y , . = y . + 4^  [y*ii  + 4yf  + 
n+1  n-1  3 L n+1  J n 

(19) 

where: 

xn  ' i(V  V t-> 

(20) 

“ 8(V  V Ca) 

(21) 

Consideration  of  the  Runge-Kutta  and  the  Milne  algorithms  shows  that 
K'tlr.e'et  method  requires  only  two  evaluations  of  f and  g per  time  step, 
•whereas.  Runge-Kutta  requires  four.  This  makes  Milne's  method  more  ef- 
ficient. 

Steady  flow  between  a source/sink  pair.  Consider  the  steady  flow  of- 
fluid  between  a line  source  and  a line  sink  of  equal  strength  * Q, 
separated  by  a distance  ■*  a.  Let  the  source/sink  pair  be  located  in  a 
homogi-naous,  isotropic  medium,  of  infinite  extent,  and  saturated  thickness 
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b.  Suppose  that  at  time  zero  the  concentration  of  solute  at  the  source  is 
changed  from  zero  to  C0. 


The  well-known  tine  dependent  solution  for  the  concentration  of  fluid 
recovered  at  the  sink  depends  on  the  travel  time  of  fluid  particles  (e.g. 
Charbeneau  and  Street,  1979).  Travel  time  t is  expressed  as  a function  of 
0,  the  direction  of  travel  of  a particle  as  it  issues  from  the  source.  If 
the  angle  8 is  measured  from  a line  between  the  source  and  sink,  then: 

c - -^-20  [1  - 9 cote]  e > o (22) 


where  <j>  is  the  effective  porosity.  For  0 = 0,  the  minimal  travel  time  t^, 
is: 


irb<fra2 
m “ 3Q  * 

At  time  t > tm  the  relative  concentration  of  recovered  fluid  is: 


(23) 


(24) 


where  the  function  0(t)  is  defined  by  (22). 


The  analytical  solution  to  the  source/sink  problem  was  compared  with  a 
numerical  solution  obtained  via  the  method  of  characteristics  and  the 
Runge-Kutta  algorithm.  Variables  were  assigned  the  values:  a » 500  m,  Q * 
10000  m3/days,  b 3 50  a,  <}>  » 0.2.  The  numerical  calculation  began  with  At 
■ 0.05  day;  subsequent  values  of  At  were  selected  so  as  to  allow  fluid 
particles  to  travel  about  10  meters  per  time  step.  Both  analytical  and 
numerical  results  are  plotted  in  Figure  1;  note  that  the  two  solutions 
practically  coincide. 


Since  the  concentration  of  fluid  recovered  at  the  3ink  varies  with 
time,  the  transport  may  be  viewed  as  a mixing  process.  In  fact,  this 
mixing  or  "apparent  dispersion"  is  obviously  just  the  result  of  nonuniform 
flow.  The  source/sink  example  leads  to  the  conclusion  that  unrepresented 
nonuniforra  flow  may  result  in  considerable  unexplained  "dispersion.” 


Transport  in  heterogeneous  media.  Consider  the  transport  of  a conser- 
vative solute  in  a two  dimensional  steady  flow  field.  The  flow  domain  is 
assumed  to  be  rectangular,  L by  H,  having  the  distributions  of  transmis- 
sivity and  effective  porosity  as  shown  in  Figures  2 and  3. 

Solution  of  the  transport  problem  begins  with  the  determination  of 
flow  field,  which  in  turn  begins  with  finding  the  hydraulic  potential. 

Using  the  linearized  Boussinesq  equation,  Daly  and  Morel-Seytoux  (1981 ) 
determined  an  analytic  solution  to  this  problem  subject  to  the  boundary 
conditions  on  the  potential  h;  - 

b(0,y)  - 0 h(L,y)  -■  X 
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y Distance 


17  (x»0)  = 17  (x>H)  = 0 

Their  solution  is: 


(25) 


h<*,y)  “ fjj  l A(n,0)  sin  ^ 

n=l 

+ 1 l A(n,m)  sin  — - cos  ^p-  (26) 

m=l  n=l 

where  the  Fourier  coefficients  A,  dependent  on  taedium  heterogeneity,  are 
found  by  the  application  of  an  integral  transform  method. 

The  specific  discharge  (and  average  linear  velocity)  associated  with 
the  potential  distribution  of  (26)  is  obtained  from  Darcy's  law.  Equation 
(26)  is  easily  differentiated  to  yield  the  hydraulic  gradient.  For  the 
problem  presented  here:  L * 2500  m,  H = 1800  m,  A 3 10  m,  and  saturated 

thickness  is  assumed  constant  and  equal  to  50  m.  Figure  4 is  a vector  dia- 
gram of  specific  discharge. 

The  trajectories  of  fluid  particles  (located  initially  along  the  right 
hand  vertical  edge  of  Figure  4)  were  calculated  by  the  Runge-Kutta,  predic- 
tor-corrector method.  A constant  time  interval  of  100  days  was  used.  The 
calculated  trajectories  define  the  streamlines  shown  in  Figure  5;  triangles 
are  used  to  locate  particles  at  100-day  intervals. 

The  movement  of  e sharp  concentration  front  through  the  medium  is 
shown  in  Figure  6.  It  is  assumed  that  at  time  zero  the  concentration  of 
fluid  along  the  right  hand  boundary  was  instantaneously  changed  from  zero 
to  CQ.  In  the  figure,  the  front  is  plotted  at  200-day  intervals.  The 
movement  of  any  particular  point  on  the  front  is  found  by  following  that 
point  along  its  associated  streamline. 

Consider  the  fluid  which  exits  the  porou;;  medium  at  the  left  hand  side 

of  Figure  4.  The  average  concentration  of  that  fluid  can  be  determined  by 

calculating  the  time  of  breakthrough  of  many  individual  stream  tubes.  The 
ratio  of  the  outflow  produced  by  the  tubes  which  have  broken  through  to  the 

total  outflow  can  be  obtained  at  any  time.  That  ratio  gives  the  relative 

concentration  of  the  fluid  flowing  out  of  the  medium.  Using  many  more 
stream  tubes  than  are  shown  in  Figure  5,  the  calculation  procedure  was  per- 
formed. The  result  Is  plotted  as  the  solid  line  in  Figure  7 (the  "observ- 
ed" breakthrough  curve). 

Using  dispersion  to  account  for  the  shape  of  the  breakthrough  curve. 
Suppose  that  the  existing  heterogeneity  of  the  preceding  problem  is 
unknown.  Suppose  also  that  an  experiment  is  conducted  to  determine  a 
breakthrough  curve;  the  result  is  the  "observed"  breakthough  curve  of 
Figure  7. 


82 


x Distance  Cm) 


Figure  4 


Specific  discharge;  flow  pattern  for  two  dimensional  steady 
state  Clow. 
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Figure  6.  Movement  of  a solute  front  in  the  flo,;  field  of  Figure  4 
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Figure  7,  Breakthrough  curves  for  the 


moving  solute  front  of  Figure  6 


If  the  porous  medium  were  assumed  homogeneous,  all  streamlines  in 
Figure  5 would  be  parallel  to  the  x axis.  Flow  would  be  steady  and  uniform 
and  the  problem  could  be  considered  one  dimensional.  A solution  to  the  one 
dimensional  advection-dispersion  equation  for  steady  flow  has  been  obtained 
by  Ogata  and  Banks  (1961).  The  differential  equation  is: 


JJC 

3t 


- D 


, 32C 
"Sx2 


3C  n 
+ u — = 0 
3x 


u > 0 


(27) 


subject  to: 


C(0,t)  = C0;  C(«,t)  = 0 

and  the  Initial  condition: 


(28) 


C(x,0)  * 0 


(29) 


where  u =*  average  linear  velocity;  D'  = D/<{>.  After  a change  of  origin,  the 
solution  for  the  breakthrough  curve  is:  at  x = 0: 


C 1 = L-ut  , 

— - 2 erf c + e 

o 


uL 

F* 


erfc 


L+ut 


(30) 


2500  meters  and 


(31) 


2/D’t  2/D’T 

For  the  problem  considered  here  L 

u *»  — =*  0.6.745  m/day 
9 

Several  estimates  of  the  coefficient  D'  can  be  made.  The  resulting  break- 
through curves  are  plotted  in  Figure  7 for  the  estimates  D'  = 50  and  100 
m2/day.  Note  that  these  two  curves  give  an  approximate  fit  to  the  observed 
breakthrough  curve. 

In  a series  of  experiments  dealing  with  one  dimensional  dispersion, 
Harleman  et  al.  (1963)  correlated  dispersion  coefficient  with  flow  and 
media  properties.  A variety  of  unconsolidated  materials  were  used.  The 
flow  and  transport  problem  were  such  that  the  analysis  of  Ogata  and  Banks 
(1961)  could  be  applied.  Determining  the  breakthrough  curve  and  the  aver- 
age linear  velocity  gave  Harleman  et  al.  the  ability  to  estimate  D*  from 
Equation  (30).  Their  correlation  formula  predicts  for  sand  grains: 


- - 0.90  (R,  )1,2 

v dso' 

where:  v is  the  kinematic  viscosity  [L2/T], 

Md50 

R,  “ 

d50  v 

and  d5Q  is  the  50%  grain  size  of  the  porous  material. 


(32) 


(33) 


If  the  shape  of  the  observed  breakthrough  curve  of  Figure  7 is  assumed 
to  be  the  result  of  dispersion,  (32)  can  be  u3ed  to  estimate  d50  for  the 
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2/j  a sz  &Q  'M  nr  for  D1  85  50  ffi  /do. y* 
porous  material.  «lth  O'  - 100  »2/d.y,  d50  49.32  », 

d50  = 27.68  m. 

conclusions.  The  above  results  J”"“‘““yca2ot  S’properly'Ic- 

r«  ^ ««*»*- «. 

by  the  effects  of  heterogeneity, 
suits  through  project  P.387.05.0056. 
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